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Abstract 

We compare the trace anomaly, strangeness and baryon number fluctuations calculated in lattice 
QCD with expectations based on hadron resonance gas model. We find that there is a significant 
discrepancy between the hadron resonance gas and the lattice data. This discrepancy is largely 
reduced if the hadron spectrum is modified to take into account the larger values of the quark 
mass used in lattice calculations as well as the finite lattice spacing errors. We also give a simple 
parametrization of QCD equation of state, which combines hadron resonance gas at low tem- 
peratures with lattice QCD at high temperatures. We compare this parametrization with other 
parametrizations of the equation of state used in hydrodynamical models and discuss differences 
in hydrodynamic flow for different equations of state. 
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I. INTRODUCTION 



Equation of state of hot strongly interacting matter play important role in cosmology [l|, |2| 
and hydrodynamic description of heavy ion collisions [3]. In many cases hydro dynamical 
models which try to describe the collective flow in heavy ion collisions used equation of state 
(EoS) with first order phase transition, although lattice QCD shows that the transition to 
the deconfined phase is only a crossover j^]. It is not obvious to what extent the collective 
flow is sensitive to details of the equation of state (EoS). It turns out, however, that in 
ideal fluid dynamics the anisotropy of the proton flow is particularly sensitive to the QCD 
equation of state ^ and using lattice inspired EoS with crossover transition overpredicts 
proton elliptic flow. 

Attempts to calculate EoS on the lattice have been made over the last 20 years (see Refs. 
@, 0] for reviews). One of the difficulties in calculating EoS on the lattice is its sensitivity 
to high momentum modes and thus to the effects of finite lattice spacing. This problem is 
particularly severe in the high temperature limit. Therefore in recent years calculations have 
been done using improved staggered fermions with higher order discretization of the lattice 
Dirac operators [sl, ^ which largely reduce the cutoff dependence in the high temperature 
region. However, there is another source of discretization effects if staggered quarks are used. 
The staggered fermion formulation does not preserve the flavor symmetry of continuum 
QCD. Because of this the spectrum of low lying hadronic states is distorted and this may 
effect thermodynamic quantities in the low temperature region. 

Hadron resonance gas (HRG) turned out to be very successful in describing particle abun- 
dances produced in heavy ion collisions [13] ■ It was also used to estimate QCD transport 
coefficients ll[ as well as chemical equilibration rates [lH close to the transition temper- 
ature. Thermodynamic quantities calculated in lattice QCD with rather large quark mass 
agree well with the HRG model if the masses of the hadrons in the model are tuned ap- 



propriately to match the large quark mass used in lattice calculations [13|. Furthermore, 
the ratio of certain charge susceptibilities are not very sensitive to the details of the hadron 
spectrum and the lattice calculations of these ratios show a reasonably good agreement with 
HRG model predictions at low temperatures [Tl- 16|. The purpose of this paper is to con- 
front the results of recent lattice calculations performed with light quark masses with the 
prediction of the HRG model and clarify its range of applicability. As we will see the HRG 
model describes thermodynamic quantities quite well up to unexpectedly high temperatures. 
Therefore lattice EoS can be combined with HRG EoS at low temperatures to get rid of 
large discretization effects. Such an EoS is also useful for hydrodynamic modeling, and we 
construct a parametrization for an EoS interpolating between HRG at low temperatures and 
lattice QCD at high temperatures. The rest of the paper is organized as follows. In section 
II we discuss the hadron resonance gas model and the effects of finite lattice spacing on 
hadron masses. Section HI deals with the comparison of fiuctuations of baryon number and 
strangeness with the prediction of HRG model. In section IV we compare the HRG with 
the lattice results on trace anomaly and construct a parametrization of equation of state 
which is easy to use in hydrodynamic simulations. In this section we also give a detailed 
comparison with other parametrizations of EoS in the literature. In section V we discuss 
hydrodynamic fiow for different EoSs discussed in the previous section, and highlight the 
differences. Finally section VI contains our conclusions. In appendices we discuss some 
technical aspects of the calculations, in particular the fits of the lattice data as well as the 
hydrodynamic fiow for partial chemical equilibrium. 
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II. HADRON RESONANCE GAS AND LATTICE QCD 



At sufficiently low temperature thermodynamics of strongly interacting matter at zero 
baryon number density is dominated by pions. The interaction between pions is suppressed 
and chiral perturbation theory can be used to estimate the pion contribution to thermody- 
namic potential [13] ■ In fact, for temperatures T < 150 MeV, the energy density of pions 
calculated in 3-loop chiral perturbation theory differs only by less than 15% from the ideal 
gas value [13] • As temperature increases heavier hadrons start to contribute to thermody- 
namics. For temperatures T > 150MeV heavy states dominate the energy density. However, 
the densities of heavy particles are still small, rii ~ exp(— Mj/T), and their mutual interac- 
tions, being proportional to rijnfc ~ exp( — (Mj + Mk)/T), are suppressed. Therefore one can 
use the virial expansion to evaluate the effect of interactions |l8|. I n the low temperature 
limit the virial expansion reduces to chiral perturbation theory [l7|. The virial expansion 
together with experimental information on the scattering phase shifts was used by Prakash 
and Venugopalan to study thermodynamics of low temperature hadronic matter [19!]. Their 
analysis showed that there is an interplay between attractive interactions (characterized 
by positive phase shifts) and repulsive interactions (characterized by negative phase shifts) 
such that their net effect can be well approximated by inclusion of free resonances : p, 
K*, A(1234) etc. Thus the interacting gas of hadrons can be fairly well approximated by 
a non-interacting gas of resonances corroborating earlier ideas of the statistical bootstrap 
model j20|. To summarize, the partition function of strongly interacting matter at low tem- 
perature can be well approximated by the partition function of non-interacting hadrons and 
resonances 



P 



j£ mesons 

+^ 5^ lnZ^^(T,r,/ixO, (2.1) 

i£ baryons 



where 



In = / dkP ln(l T ^.e"-/^) , (2.2) 

with energies Ei = \/k^ ^ mf, degeneracy factors di and fugacities 

z, = exp|^(^X>xO/2^j • (2.3) 

Here we consider all possible conserved charges X"-, including the baryon number B, electric 
charge Q, strangeness S etc. We do not include any repulsive interactions in the form of 
excluded volume corrections or repulsive mean field |2l[ ]. 

The assumption that thermodynamics in the low temperature region is well described 
by a gas of non-interacting hadrons and resonances is important for practical applications 
of hydrodynamic models. At the end of the hydrodynamical evolution, the fluid is usually 
converted into particles using the Cooper- Frye procedure [22j . This procedure conserves en- 
ergy, momentum and charge without any specific considerations if, and only if, the equation 
of state of the fluid is the same than the equation of state of free particles [23|. There- 
fore it is important to confront the predictions of the hadron resonance gas for different 
thermodynamical quantities with the available lattice data. 
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The most extensive lattice calculations of the equation of state have been performed 
using two versions of improved staggered fermions: the so-called asqtad and p4 formulations 
[sl, M, US- Calculations have been performed using lattices with temporal extent A^,- = 4 
and 6 (sl, 0] and more recently with temporal extent Nt- = 8 [24]. These correspond to 
lattice spacings a = 1/(4T), 1/{6T) and 1/(8T) respectively. The strange quark mass 
was close to its physical value, while the light {u and d) quark masses were one tenth of the 
strange quark mass, rrig = Clm^, corresponding to the pion mass in the range 220 — 260 
MeV. The lattice calculations of the equation of state have been compared to the prediction 
of the hadron resonance gas jo], H^l ■ It turned out that the lattice results fall considerably 
below the HRG prediction. One obvious reason for this discrepancy is the fact that the 
quark mass used in lattice calculations is about a factor two larger then the physical one. 
However, this fact alone is unlikely to explain the whole discrepancy as the contribution 
of the pseudo-scalar mesons to the energy density is small and quark mass dependence of 
other hadron masses in this small quark mass region is relatively weak. On the other hand 
the lattice spacing dependence of the hadron masses may play an important role. Since the 
lattice calculations of the EoS are performed at fixed temporal extent A^^, the temperature 
is varied by changing the lattice spacing T = l/{Nra). As the temperature is decreased the 
lattice spacing gets larger and the cutoff effects on the hadron masses increase, i.e. the size 
of cutoff effects on the hadron masses is a function of the temperature. 

The hadron masses for asqtad improved staggered fermions have been studied in Refs. 



25H28| for several lattice spacings a ~ 0.06, 0.09, 0.125, 0.18 and 0.25 fm. For all of these 



lattice spacings there are significant deviations in the hadron masses from the experimental 
values. Of course, after the proper continuum extrapolations all the hadron masses agree 



with the experiment j25|]. The large cutoff dependence of the hadron masses in the staggered 
formulation is due to large 0{asa'^) discretization errors which also break the flavor symme- 



try. This is not the case for the improved Wilson fermion formulation [29|, |30| , where cutoff 
dependence of the hadron masses is below 5% already at lattice spacing < 0.2fm. In the 
following subsections we discuss the cutoff dependence of the pseudo-scalar meson, vector 
meson and baryon masses separately. 



A. Pseudo-scalar mesons in staggered formulation 

The staggered formulation of lattice QCD describes four degenerate quark flavors in the 
continuum limit. To obtain the physical number of flavors, i.e. one relatively heavy strange 
quark and two light quarks, the so-called rooting procedure is used. The rooting procedure 
amounts to replacing the fermion determinant in the path integral expression of the partition 
function with its fourth root. In the continuum limit this procedure is justifled j3l|^. At flnite 
lattice spacing, however, the four flavors are not degenerate, but there are flavor changing 
discretization effects of order O^asO^)- As the result of this the sixteen pseudo-scalar mesons 
of the 4 flavor theory have unequal masses, and only one of them has a vanishing mass in 
the limit of the zero quark mass, ruq — )■ 0. The sixteen pseudo-scalar mesons can be grouped 
into eight multiplets js^l- The multiplets are characterized by the different masses mpg. 
and degeneracies d^^, i = 0,1, ..7. The flrst multiplet contains one Goldstone pseudo-scalar 



^ The justification of tlie rooting procedure at finite lattice spacing is still subject of debate see e.g. Ref. 
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TABLE I: The parameters of Eq. (|2.5p describing the quadratic pseudo-scalar meson sphttings. 
Also shown are the flavor matrices and the degeneracy factors d*. 



meson, i.e. '^ps^^ ~ and c/p^ = 1. The flavor generators of the 4 flavor theory can be 



chosen to be the product of the Dirac matrices [34, 135| . The masses of other pseudo-scalar 
mesons are given by 

"^Js, = "^Jso +^"^Js,- (2-4) 

The quadratic splittings ^'rv?ps^ are independent of the input quark mass to a very good 
approximation and are proportional to (agO^) for small lattice spacings, and in general, 
increase with increasing index %. The correspondence between the index i and the flavor 
matrix as well as the degeneracies rfpg for non-Goldstone pseudo-scalar mesons are given in 
Table [11 

The quadratic splittings have been calculated in numerical simulations with asqtad action 



in Refs. |25|, |27[. For lattice spacings 0.09fm < a < 0.25fm the data can be well parametrized 
by the form 

-I ■ Hs. = "ir^-^^ - = (-/rif- (2-5) 

Here and in what follows we use the scale parameter ri extracted from the static quark 
potential V{r) to convert from lattice units to physical units. The scale parameter ri is 
defined as 



2 



dV{r) 
dr 



I r=ri 



1.0. (2.6) 



We use the value ri = 0.318 fm determined from bottomonium splitting [27|. The values 
of the parameters Op^, 6pg, Cpg, /Spg are given in Table [B In Figure [J we show the quadratic 



pseudo-scalar meson splittings calculated by the MILC collaboration [25|, |27| as function of 
the lattice spacing and compared with the parametrization given by Eq. ( 12. 5p . As one can 
see from the Figure [H this parametrization gives good description of the data. In Figure 
[1] we also show the pseudo-scalar meson splitting for stout action used by the Budapest- 
Wuppertal group [36 1. The splittings are significantly reduced compared to the calculations 
with asqtad action. To take into account the fiavor symmetry breaking in the pseudo-scalar 
meson sector, i.e. the fact that pion masses are non-degenerate the contributions of pions 
and kaons to the pressure is calculated as 

P^'^/^' = ^T7^ E^;sln2^^(mp,), (2.7) 

i=0 
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FIG. 1: The quadratic splittings of non-Goldstone pseudo-scalar mesons in the seven different 
multiplets calculated with asqtad action 25|] at different lattice spacings. The lines show the 
parametrization given by Eq. (j2.5p . The open symbols refer to the lattice data obtained with the 
stout action [36.] . 



where mpg., i = 1 — 7 is calculated according to Eq. (12. 4 p and mps^ is equal to the pion or 
kaon mass used in the actual lattice calculations. Figure [T] shows that the splitting between 
different pseudo-scalar mesons is quite large for lattice spacings used in calculations of the 
equations of state (a = 0.12 — 0.18)fm. Even in the calculations with the stout action the 
mass of the heaviest pion that enters Eq. (12. 7p is about SOOMeV if A^,- = 8 lattices are used. 
As the result the contribution of the pseudo-scalar mesons to thermodynamic quantities at 
finite lattice spacings is smaller than in the continuum. 



B. Vector meson masses 



In lattice calculations hadron masses are functions of the input quark mass and the lattice 
spacing. The quark mass dependence of the hadron masses is usually studied in terms of the 
lightest (Goldstone) pion mass. Since in lattice calculations quark masses are larger than 
the physical ones, and there is no full flavored chiral symmetry at finite lattice spacing a 
combined extrapolation in the pion mass and lattice spacing is needed for all hadron masses. 
We fitted the lattice spacing and pion mass dependence of the vector meson masses with 
the following simple formula 

n ■ ,n(a, ™j = n™*' + + j^. . = (a/nf, (2.8) 

1 + a2X 1 + 02X 

where m^^^* is the physical value of the meson mass from the particle data book. It turns 
out that this formula describes the lattice spacing dependence of the vector meson masses 
in the interval 0.06fm < a < O.lSfm for (rim^)^ < 0.8. For strange vector mesons we set 
02 = as there is no apparent lattice spacing dependence of the slope in their quark mass 
dependence. The numerical values of ai, 02, &i and 62 are given in the Appendix A, where 
also the lattice data on the vector meson masses are shown. 
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C. Baryon masses 



The nucleon and the baryon masses have been calculated by the MILC collaboration 



with asqtad action at five different lattice spacings, a = 0.06,0.09,0.12 and O.lSfm |25|, [27 



28l . I37j and several quark masses. We performed a simultaneous fit of their quark mass 
(pion mass) and the lattice spacing dependence using the Ansatz given by Eq. (12.81) . which 
works well for (rimTr)^ < 0.8. The parameters of this fit together with the lattice data on 
the nucleon and Q baryon masses are presented in the Appendix A. In calculations with 
asqtad action the value of the strange quark mass was slightly larger than its physical value 
for each lattice spacing considered. This is due to the fact that the strange quark mass was 
fixed considering the ratio of the meson mass to the mass of the unmixed rjss pseudo-scalar 
meson instead fixing the kaon mass. This gives difference O(a^) in the value of the strange 
quark mass. To take this into account the lattice spacing dependence of the strange quark 
mass was parametrized as ms/m^'^^*(a) = 1 + 1.02(a/ri)^. Then to estimate the Q baryon 
mass we used the following formula 

rimn(a,m^) = rim^^^ + ai{rimn)'^ — ai{rim^'^''^^)'^ + bix+{m^^^ — m'^^^'^)-1.02x, x = {a/rif'. 

(2.9) 

Here the last term accounts for the small deviation of the strange quark mass from its 
physical value. For other baryons (A, A, S, and H) no such detailed lattice calculations 
are available. Therefore to estimate the cutoff dependence of the A mass we use the same 
formula as for nucleon. While to estimate the cutoff dependence of the masses of the strange 
baryons we used the following formulas 

3 1 + a2X 1 + b2X 1 + a2X \ml / 

r, . m^ia, mj = m, + 5^^^ + + [-^j (2.11) 

3 l + a2X I + O2X 1 + 020; \mT J 

X = (a/riY. 

Here again we have taken into account that the strange quark mass in simulations with 
asqtad was slightly larger than the physical value. In the Appendix we give the comparison 
of the baryon masses calculated using the above formulas with available lattice results. It 
turns out that our parametrization of the baryon masses works reasonably well. We will 
use these parametrizations when calculating different quantities in the HRG model in the 
following sections. 



III. FLUCTUATIONS OF CONSERVED CHARGES 

Derivatives of the pressure with respect to chemical potential of conserved charges, e.g. 
baryon number (B), electric charge Q and strangeness 5* can be easily calculated in lattice 
QCD 

Xn - J- Q-f, Ux=0, A - ^. [6.1) 
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These are related to quadratic and higher order fluctuations of conserved charges X2 — 
{X^)/{VT^), xf = {{Nj.) - 3{N]^)^)/{VT^) etc.^ Contrary to the pressure itself the evalu- 
ation of these derivatives does not involve zero temperature subtraction and integration in 
the temperature variable starting from some low temperature value. Therefore it is easy to 
compare them to the prediction of the HRG model. Different fluctuations up to the sixth 



order have been calculated for p4 and asqtad action pj|, 138|, 139| on A^^ = 4 and A^^ = 6 
lattices. Quadratic strangeness fluctuations have been also calculated on A^,- = 8 lattices 
for the p4 and asqtad actions by the HotQCD collaboration [2^]. While there are extensive 
calculations with p4 action at flnite temperature, the zero temperature hadron spectrum was 
not studied in detail for p4 action. Therefore our analysis of fluctuations and comparison 
with HRG model will mostly rely on results obtained with asqtad action. 

Let us start our discussion with baryon number fluctuations. Baryon number fluctuations 



for asqtad action have been calculated in Ref. [38| for nig = 0.2ms on = 6 lattices. 
Equations (I2.8p - fl2.12p describe the quark mass and lattice spacing dependence of ground 
state baryon masses calculated on lattice with asqtad action. Nothing is known about the 
lattice spacing dependence of the excited baryon masses which play very important role in 
baryon number fluctuations in the temperature range of interest. We assume that all the 
excited baryons up to certain mass threshold mf^f- have the same quark mass and lattice 
spacing dependence, while the baryon masses above that threshold are not modifled by flnite 
lattice spacing. The mass threshold m^^^ is an additional parameter of our model. In Fig. 
[2] we show the lattice data for baryon number fluctuations for asqtad action compared with 
hadron resonance gas model with physical value of the baryon masses including all the states 
up to 2.5GeV. The lattice data fall considerably below the HRG prediction. The baryon 
number fluctuations have also been calculated in a HRG model, where all the baryon masses 
up to the Tn^^-f- = 1. 8GeV and 

~ 2.5GeV have been modifled according to Eqs.(|2^- 
fl2.12p . The corresponding results are shown as dashed lines. As one can see, the HRG 
overshoots the lattice data with = 1.8GeV, while with = 2.5GeV it undershoots 
the lattice data. However, the agreement between lattice data and HRG is greatly improved. 
For completeness we also show the lattice data for p4 action calculated for rrig = 0.1ms (l6| . 

Strangeness fluctuations have also been calculated on the lattice using asqtad and p4 



action [16|, |24, |38|, |39[ on A^,- = 6 and A^r = 8 lattices. We have calculated strangeness 
fluctuations in HRG model, where ground vector state meson masses have been calculated 
using Eq. (12. 8p . while baryon masses have been calculated using Eqs. fl2.8p - fl2.12p . The con- 
tribution of kaons has been treated in the way discussed in section Hi At i.e. for each physical 
kaon averaging over sixteen staggered flavors have been performed and the mass splitting 
has been calculated using Eq. (12. 5p . This turns out to be important for the description of 
strangeness fluctuation for temperatures T < 165MeV, where the contribution of kaons is 
quite signiflcant. In Fig. [3] we show the prediction of the HRG model with modifled hadron 
masses compared to the lattice data for asqtad action for A^^ = 8. In the flgure we show the 
prediction of the HRG model with cutoffs = 1.8GeV and = 2.5GeV for the baryon 
masses. Here the effect of using different cutoffs for the modiflcation of the baryon masses is 
signiflcantly smaller as meson contribution to strangeness fluctuations is large. In the flgure 
we also show the prediction of the HRG model with physical hadron masses as well as the 
lattice results for the p4 action. As one can see the lattice data fall significantly below the 



Here we consider the case of zero chemical potential, so (Nx) = 0. 
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FIG. 2: Baryon number fluctuations calculated with asqtad action on the N-r = 6 lattices compared 
with HRG model with physical value of the baryon masses (solid line) and with HRG model with 
baryon masses calculated according Eqs. (|2.8p - (|2.12p . m^^^ = l.SGeV and = 2.5GeV (dashed 
lines). Also shown are the lattice results for the p4 action. 

predictions of HRG model with physical quark masses, while the HRG model with modified 
hadron masses gives quite a good description of the lattice data. For completeness we also 
show the result for HRG model with modified hadron masses for A'^ = 12. As one can see 
from Figs. |2] and |31 HRG model can describe baryon number and strangeness fiuctuations 
reasonably well up to temperatures as high as Tc. 

IV. QCD EQUATION OF STATE 

A. The trace anomaly and parametrization of the equation of state 

Available lattice data provide an EoS which is not easy to use in hydrodynamic models 
because different thermodynamic quantities are suppressed in the low temperature region 
due to discretization errors. In the previous section we have seen this for baryon number 
and strangeness fiuctuations. 

In lattice QCD the calculation of the pressure, energy density and entropy density usually 
proceeds through the calculation of the trace anomaly 6(T) = e(T) — 3p(T). Using the 
thermodynamic identity the pressure difference at temperatures T and Tiow can be expressed 
as the integral of the trace anomaly 



By choosing the lower integration limit sufficiently small, p(Tiow) can be neglected due to 
the exponential suppression. Then the energy density e(T) = 0(T) + 3p(T) and the entropy 
density s(T) = (e+p)/T can be calculated. This procedure is known as the integral method 




(4.1) 
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FIG. 3: The strangeness fluctuations calculated on A',- = 8 lattices for asqtad and p4 actions and 
compared with the prediction of the HRG model with physical (solid line) and modified (dashed 
lines) hadron masses. The upper (lower) dashed line corresponds to = 1.8(2.5)GeV. The 

doted lines show the prediction of the HRG with modified hadron masses for A^^^ = 12. 



40[. Since the trace anomaly plays a central role in lattice determination of the equation of 
state, we will discuss it in the HRG model and its comparison with the lattice data in the 
following. As we will see this helps constructing realistic equation of state that can be used 
in hydrodynamic models. 

As mentioned before, finite temperature lattice calculations are usually performed at 
fixed temporal extent A^^ and the temperature is varied by varying the lattice spacing a, 
T = l/(iV^a). Thus, calculations at low temperatures are performed on coarse lattices, 
while the lattice spacing gets smaller as the temperature is increased. Consequently the 
trace anomaly can be accurately calculated in the high temperature region, while in the 
low temperature region it is affected by possibly large discretization effects. Therefore to 
construct realistic equation of state we could use the lattice data for the trace anomaly in 
the high temperature region, T > 250MeV, and use HRG model in the low temperature 
region, T < 180 MeV. In Fig. H] we compare the lattice results on trace anomaly obtained 
on A^^ = 8 lattices with asqtad and p4 action with the HRG model. The HRG model with 
modified masses appears to describe the lattice data quite well up to temperatures of about 
180MeV. In the intermediate temperature region 180MeV < T < 250MeV the HRG model 
is no longer reliable, whereas discretization effects in lattice calculations could be large. The 
later can be seen by comparing the lattice data obtained on A^^ = 6 and A^^ = 8 lattices 
with p4 and asqtad action. Therefore we constrain the trace anomaly in the intermediate 
region only by the value of the entropy density at high temperatures. 

In pure gauge theory, where continuum extrapolation has been performed, the entropy 
density falls below the ideal gas limit only by 15% at temperatures of about 4Tc jioj [T^ is the 
transition temperature). In QCD the entropy density calculated on A^^ = 6 and 8 lattices is 
(5 — 10)% below the ideal gas limit j4l[ at the highest available temperature. Furthermore, 
resummed perturbative calculations describe quite well the entropy density in the high 
temperature region both in pure gauge theory j42[ and in QCD [41]. These calculations 
indicate that deviation from the ideal gas limit is about (5 — 10)%. The fluctuations of 
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FIG. 4: The trace anomaly calculated in lattice QCD compared with the HRG model with physical 
hadron masses (solid line) and modified hadron masses (dashed lines). The upper (lower) dashed 
line corresponds to = 1.8(2.5)GeV. 



quark number are also close to ideal gas limit and well described by resummed perturbative 



calculations j4l[. Here the deviation from the ideal gas limit is less than 5%. Therefore we 
will use the guidance from existing lattice QCD calculations and require that the entropy 
density is below the ideal gas limit by either 5% or 10%, when parametrizing the trace 
anomaly. 

At high temperature the_trace anomaly can be well parametrized by the inverse poly- 
nomial form (see e.g. Ref. [2^). Therefore we will use the following Ansatz for the high 
temperature region 

(e - 3P)/T^ = c/a/T^ + di/T^ + Ci/T^^ + Ca/T'^^ (4.2) 

This form does not have the right asymptotic behavior in the high temperature region, where 
we expect (e — 3P)/T^ ~ Q^^iT) 1/ \n^{T/ Aqcd), but works well in the temperature range 
of interest. Furthermore, it is flexible enough to do the matching to the HRG result in low 
temperature region. We match to the HRG model at temperature Tq by requiring that the 
trace anomaly as well as its first and second derivatives are continuous. The parametrization 
of the trace anomaly and thus QCD equation of state obtained using these requirements are 
labeled by s95p— vl, s95?2— vl and s90/— vl. The labels "s95" and "s90" refer to the fraction 
of the ideal entropy density reached at T = SOOMeV ( 95% and 90% respectively), whereas 
the labels p, n and / refer to a specific treatment of the peak of the trace anomaly or its 
matching to the HRG. The detailed procedure of performing the fit to the lattice data and 
matching to the HRG model is described in Appendix B, where also the labeling scheme is 
explained in more detail. The values of the parameters Tq, d2, d^, ci, C2, rii and n2 in each 
parametrization are given in Table [Tll The HRG result for the trace anomaly can also be 
parametrized by the simple form 

e — 3P 

-j;^ = aiT + a^T"" + a^T^ + a^T^\ (4.3) 
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d2 (GeV2) 


(GeV^) 


ci (GeV"i) 


C2 (GeV"2) 


ni n2 


To (MeV) 


s95p 


0.2660 


2.403-10~3 


-2.809-10-^ 


6.073-10-23 


10 30 


183.8 


s95n 


0.2654 


6.563-10-3 


-4.370-10-5 


5.774-10-6 


8 9 


171.8 


s90/ 


0.2495 


1.355-10-2 


-3.237-10-3 


1.439-10-^^ 


5 18 


170.0 



TABLE II: The values of the parameters for different fits of the trace anomaly. 



with ai = 4.654 GeV-\ = -879 GeV-^, = 8081 GeV-^ = -7039000 GeV-^" (see 
Appendix B for details). 

The lattice data for trace anomaly compared to the parametrization given by Eqs. fl4.2p 
and (14. 3 p is shown in Fig. |5l We show three parametrizations in the figure corresponding to 
a entropy density at T = 800MeV which is below the ideal gas limit by 5% (s95p— vl and 
s95n— vl, the solid and dotted lines, respectively) and 10% (s90/— vl, dashed line). All the 
parametrizations describe the lattice data for T > 250MeV, while in the low temperature 
region, T < 170MeV, they are significantly above the lattice results. On the other hand, our 
parametrizations of the trace anomaly are below the lattice data in the peak region. This 
comes from the imposed constraint on the entropy density at T = 800MeV. If we would use a 
parametrization which goes through the N^- = 8, p4 lattice data and matches the resonance 
gas model at some temperature near 190MeV, the entropy density would overshoot the ideal 
gas limit already at temperatures of about 600 — 700MeV. Such a behavior would contradict 
the available lattice and weak coupling results. 

The difference in the s95n— vl and s95p— vl parametrizations is in the treatment of the 
peak region. When we do the fit only on the lattice data above T = 250MeV (s95n— vl, 
dotted line), the peak value is clearly below the present data. To explore the sensitivity of the 
EoS to the height of the peak, we also did the fit using one additional point at T = 205MeV 
close to the present data, and the same entropy constrain (s95p— vl, solid line). This forces 
the trace anomaly almost to reach the data at the peak maintaining a reasonable fit to the 
data at high temperatures. 

The EoSs, obtained by integrating the parametrizations given in Eqs. (14. 2 p and (14.30 over 
temperature as shown in Eqs. (14.10 . are shown in Fig. [6l The clearest difference between 
our different parametrizations is trivial: The different behavior at high temperatures is due 
to the different entropy constraint at T = 800MeV, which is of course manifested in energy 
density and pressure too. On the other hand, the different height of the peak of the trace 
anomaly causes only a tiny difference in pressure and energy density around T = 200 MeV. 
This difference is manifested much clearer in the speed of sound: The large peak and larger 
switching temperature from HRG to lattice causes much more rapid change in the speed 
of sound in s95p— vl parametrization than in the two others. One could claim that the 
differences in the speed of sound are due to the different matching temperatures Tq, but we 
want to remind the reader that we treated Tq as a fitting parameter in our model, and any 
changes in Tq would make the fit to the lattice data worse (see details in Appendix [B]). In 



Figure [6] we show the speed of sound for EoS Q from Refs. [43|, |4^: An equation of state 
with a first order phase transition at = 170MeV. Below EoS Q coincides with HRG 
EoS, while above this temperature it is given by bag equation of state with three massless 
flavors and bag constant B = (0.2447Ge^)^. Nevertheless, the most striking feature of the 
speed of sound in the proposed parametrization of the EoS is that there is no softening 
below the hadron gas. There is no region where speed of sound would be smaller than in 
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FIG. 5: The trace anomaly calculated in lattice QCD with p4 and asqtad actions on A',- = 6 and 
8 lattices compared with the parametrization given by Eqs. (j4.2p and (j4.3p . The solid, dotted 
and dashed lines correspond to parametrizations s95p— vl, s95n— vl and s90/— vl respectively, as 
discussed in the text. 
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FIG. 6: The pressure, energy density (left panel) and speed of sound (right panel) in the equations 
of state obtained from Eqs. (j4.2p and (|4.3p . The vertical lines indicate the transition region (see 
text). In the right panel we also show the speed of sound for the HRG EoS and EoS with first 
order phase transition (thin dotted) line, the EoS Q 

hadron gas, and its minimum value is that of HRG speed of sound'^. It is quite simple to 
understand why this happens: To achieve smaller speed of sound than the speed of sound in 
hadron gas, the trace anomaly should be larger than in HRG. As one can see in Fig. HJ the 
present lattice data clearly disfavors such a scenario. In Figure [H] we indicate the transition 
region from hadronic matter to deconfined state by vertical lines. We define the transition 



Similar EoS was presented already in Refs. 



45 
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region as the temperature interval 170MeV < T < 220MeV. In view of the crossover 
nature of the finite temperature QCD transition such definition is ambiguous. The lower 
temperature limit in our definition comes form the fact that for T < 170MeV all our EoS 
agree with HRG EoS. The rapid rise in the energy density and entropy density stops roughly 
at 220MeV and starting from this temperature the variation of thermodynamic quantities is 
quite smooth 2J|. Also the Kurtosis of the baryon number becomes compatible with quark 
gas at temperatures of about 220MeV . For these reasons we have chosen it as the upper 
limit for the transition region. 



B. Comparison with other works 



The idea of using HRG in low temperatures and parametrized lattice EoS in high tem- 
peratures is by no means new. Laine and Schroder constructed QCD equation of state based 
on the effective theory approach in the high temperature region, while using the resonance 
gas equation of state in the low temperature region [i^ . In the effective theory approach the 
contribution of hard modes is treated perturbatively, while the contribution of soft modes 
is calculated using 3 dimensional lattice simulations of the effective theory called EQCD 
This parametrization has been used in recent recent viscous hydrodynamic calcula- 



tions |49| . The smooth matching to the resonance gas was done in the temperature interval 
T = 170 - ISOMeV. 

Two other parametrizations of the EoS 45, IH, HO] used lattice data of Budapest- 
Wuppertal(BW) group obtained using so-called stout staggered fermion action [5l| and 
temporal extent Nt = 4 and 6. Since the stout action does not use higher order difference 
scheme in the lattice Dirac operator, discretization effects at high temperatures are very 
large. In the ideal gas limit the pressure calculated on A^^ = 4 and 6 lattices is about twice 
the continuum value. As the consequence the lattice results of the BW group have large 
cut-off effects and overshoot the continuum ideal gas at high temperatures, see discussion 
in Ref. [52[. In an attempt to correct this problem the authors of Ref. 51 1 divided all their 
thermodynamic observables by the corresponding ideal gas value calculated on = 4 and 
6. Since cutoff effects are strongly temperature dependent this procedure overestimates cut- 
off effects in the interesting temperature region and underestimates the pressure and other 
thermodynamic observables. 

The Krakow group used the stout lattice results parametrized in Ref. 



53 



with Tc = 

167MeV and matched the speed of sound to the HRG result at that temperature |45[ 46] . The 
procedure of the Krakow group involves extracting all the other thermodynamical quantities 
from the speed of sound using the relation 



s(T) = s(To)exp 



dV 



(4.4) 



Connecting the speeds of sound in HRG and lattice leads in this procedure to a larger entropy 
density at high temperatures than given by the lattice parametrization of Ref. [53| . To make 
their EoS to fit the lattice data at T ^ IGeV, the Krakow group made the speed of sound 
smaller in the temperature region 28 < T < 118 MeV by hand [sS]. This change is below 
the expected freeze-out temperature, and thus the speed of sound which affects the actual 
hydrodjTiamical evolution is unchanged. However, since entropy density is calculated by 
integrating over the entire temperature range, entropy density is smaller than the original 
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HRG value everywhere in the range 28MeV < T < T^. More specifically, in the range 
120 < T < 160MeV, both energy and entropy densities and pressure are ~ 5% below the 
HRG value, and energy is not automatically conserved at freeze-out, see the discussion in 
section |Vl 

Heinz and Song also used the BW lattice results, but they parametrized pressure and 
temperature as function of energy density and matched them to HRG result; they used = 
172MeV [5(|. As the authors themselves note, their EoS is not exactly thermodynamically 
consistent, which leads to a violation of entropy conservation in an ideal fluid calculation. 
However, to our understanding this does not affect the qualitative studies the EoS has been 
used so far and the conclusions of those papers should be valid. 

Here a note concerning Tc in these parametrizations is in order: In Ref. j5l| thermo- 
dynamic quantities are given as function of T/Tc but the value of is not specified. At 
lattice spacings corresponding to temporal extent A^,- = 4 and N^- = 6, used in BW EoS 
calculations, has large cutoff effects and rnay deviate considerably from the continuum 
value Tc = 170(3) (4)MeV determined in Ref. jSGj (e.g. calculations of reported at Quark 
Matter 2005 on A^^ = 4 and 6 lattices gave = 189(8)MeV |55|). The best way to elimi- 
nate part of the cutoff effects in the BW EoS is to use the continuum value of the transition 
temperature Tc, which is interestingly enough appears to agree within errors with the values 
used in the phenomeno logical parametrizations discussed above. 

The parametrization of the EoS by the HotQCD Collaboration 2J] is based on a simple 
fit of the lattice results on the trace anomaly. For temperatures below 130MeV, where no 
lattice data is available, HRG values for the trace anomaly have been used and assigned 
an artificial error. The resulting parametrization is well below the HRG at temperatures 



T > 50 MeV, see discussion in Ref. [2|]. For example, at T ^ ISOMeV and T ^ 170MeV 
temperatures, pressure, energy and entropy densities are roughly 20% and 10% smaller than 
in HRG, respectively. Therefore, when one uses this parametrization, energy conservation 
at freeze-out requires additional consideration, see the discussion in section |Vl 

In the Fig. [7] we show the comparison of our parametrization for EoS with the ones 
discussed above including the trace anomaly, speed of sound, pressure and energy density as 
function of the temperature as well as the pressure as function of the energy density. From 
the figure we see that there are significant differences between different parametrizations. 
The parametrization based on BW lattice results seem to have different high temperature 
asymptotic. In the low temperature region the EoS L parametrization by Heinz and Song 
and the HotQCD parametrization differ significantly from others. It is worth noting that 
EoS L does not even try to reproduce the HRG below 130 MeV temperature. The authors 
assumed that the details of the EoS below freeze-out temperature have only a negligible effect 
on the evolution within the freeze-out surface and a faithful reproduction of the HRG is thus 
not needed. Finally, for the trace anomaly, the difference between our parametrization and 
HotQCD parametrization is limited to temperatures T < 250MeV by construction. However, 
since the pressure, energy density and entropy density is obtained using the integral method 
the differences in trace anomaly at low temperatures result in differences at all temperatures 
for these quantities. 



V. EFFECT ON HYDRODYNAMICAL FLOW 

Now we are in a position to use the lattice results on EoS in a hydrodynamical model 
and compare the results with previous approaches using first order phase transition or the 
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FIG. 7: The trace anomaly, the speed of sound, the pressure and the energy density as function 
of the temperature for different parametrizations used in hydrodynamic models. Also shown is 
the the pressure as function of the energy density. The solid black line corresponds to s95p-vl 
parametrization. 



other lattice parametrizations in the literature. We also quantify how present uncertainties 
in the EoS parametrization affect the hydrodynamical flow. For simplicity we perform our 
analysis using ideal hydrodynamics. This is also motivated by the fact that the overall value 
and temperature dependence of the QCD and HRG transport coefficients is not known and 
any attempts to parametrize them would introduce additional uncertainties in the analysis. 

As the first step we study the sensitivity of the momentum anisotropy on the EoS. This 
is the cleanest way to address the sensitivity of hydrodjTiamic flow to the EoS as additional 
coinplications due to freezeout do not enter here. The momentum anisotropy is defined 



as 



43 



_ {TxX Tyy) 



where Txx and Tyy are the diagonal transverse components of the energy-momentum tensor 
and brackets denote averaging over the entire transverse plane. In Figure [H] we show the time 
evolution of momentum anisotropy in Au+Au collisions at ^snn = 200 GeV with b = 7fm 
impact parameter for different EoSs. The left panel shows the anisotropy calculated usins 



our different parametrizations, and EoS Q from Refs. |43l. |44| . As studied in detail in Ref. |43| 



the first order phase transition causes the build up of the fiow anisotropy to stall when most 
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FIG. 8: The time evolution of the momentum anisotropy in 6 = 7 fm Au+Au colhsions using 
the EoSs developed in this paper and the old EoS with a first order phase transition (EoS Q) 
from Ref. l43l. 1441 1 (le ft p anel) and the EoS s95p— vl compared to the various lattice EoSs in the 
literature \24 , 14514471. ISOj (right panel) . In the right panel the solid black line refers to the result 
obtained with s95p— vl parametrization. The inset in the right panel shows the temperature 
evolution in the middle of the system for different EoS. The horizontal lines indicate the transition 
region. 



of the system is in the mixed phase. There is no such a structure when the transition to 
the hadronic matter is a smooth cross over, but the anisotropy increases monotonously. 
The hardness of the EoS in plasma phase is also manifested in the early behavior of the 
anisotropy. EoS Q is much harder in that region than any of the lattice EoSs studied here, 
and thus the build up of the flow anisotropy is faster. On the other hand, the mixed phase 
makes the EoS Q much softer in average during the evolution, and the final anisotropy is 
the smallest of all EoSs studied here. The speed of sound is quite similar in EoSs s95p— vl, 
s95n— vl and s90/— vl, and consequently the development of the flow anisotropy is similar. 
When the system cools, the speed of sound stays large longest for s95p— vl, but it also 
drops fastest and stays small longest for s95p— vl. These effects cancel, and the evolution of 
the anisotropy is almost identical to s95n—vl. After that argument, the largest anisotropy 
obtained using s90/— vl may look surprising, but closer inspection of Fig. E] reveals that 
s90/— vl has always either larger or equal speed of sound than s95p— vl. Thus s90/— vl is 
harder, and it should lead to a larger anisotropy than s95p— vl. 

The old wisdom has been that elliptic flow builds up quickly during the early stages of 
the evolution and is mostly build up during the plasma phase. For example, for EoS Q, 
three quarters of the final anisotropy has been built up when the center of the system 
reaches mixed phase. For lattice based EoSs this is no longer as clear: Roughly half of the 
anisotropy is built up during the transition region, after the first three fm of the evolution, 
but the hadronic contribution from below T = 170 MeV temperatures to ep is essentially 
negligible unlike in EoS Q. This difference in the evolution is partly explained by the softer 
EoS in the plasma phase — anisotropy is built up slower. Another reason is that the transition 
region in our EoSs reaches up to e ~ 3.5 GeV/fm^ energy density, whereas EoSQ reaches 
plasma phase already at e ~ 2.15 GeV/fm^ density. 

In the right panel of Fig. |8] s95p — vl is compared to other lattice EoSs in the literature. 
The differences in flow are difficult to sort out based on the speed of sound as function 
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of temperature alone, but they are easy to understand when one looks at the pressure as 
function of energy density (Fig. [7]). It can be seen that the gradient dP/de, i.e. speed of 
sound squared, is largest for the Krakow EoS, and EoS L is approximately as stiff as the 
Krakow EoS above e ~ 1 GeV/fm^ density. Thus it is not surprising that the Krakow EoS 
leads to the largest anisotropy, and that the initial build up of the anisotropy is similar for 
the Krakow EoS and EoS L. The build up of flow deviates when EoS L reaches its softest 
region which is much softer than in any other EoS, and leads to behavior reminiscent of 
EoS Q: A sudden stall in the increase of the anisotropy. Unlike in EoS Q, however, there is 
no subsequent increase in anisotropy after the soft region has been passed. Likewise, even if 
the speed of sound in the HotQCD EoS at high temperatures is equal or even larger than in 
s95p— vl or in the EoS by Laine and Schroder, the speed of sound in the transition region is 
so much smaller than in the other EoSs, that the final anisotropy is the smallest of all. As 
well, even if the speed of sound as function of temperature look different for s95p— vl and 
Laine and Schroder's EoS, as a function of energy density they are almost equal, and thus 
these EoSs lead to basically identical build up of the flow anisotropy. 

As the next step we study the sensitivity of the spectra and elliptic flow on the EoS. We 
include freeze out into the calculation described above, and use first the same freeze out 
temperature, Tfo = 125MeV, for all EoSs^. The pion and proton spectra after resonance 
decays is shown in the left panel of Fig. [9] for EoSs s95p— vl, s95n— vl, s90/— vl and 
EoS Q. As expected the new parametrization lead to flatter spectra than EoS Q, but the 
differences between the parametrizations themselves are too small to result in significant 
differences in spectra. The p^-- differential V2 of pions and protons shown in the middle 
panel are surprisingly insensitive to the EoS. The larger flow anisotropy shown in Fig. [S] 
leads to larger pt averaged ^2, but that is mostly due to flatter spectra weighting V2{pt) 
at higher pt where it is larger, than due to V2{pt) being larger. One must also remember 
that this result is obtained using the same freeze-out temperature for all the EoSs. Before 
discussing how the EoS affects elliptic flow, one has to readjust the freeze-out temperature 
to produce similar spectra. This we have done in the right panel of Fig. [9j When one uses 
Tfo = 140MeV for EoSs s95p— vl, s95n— vl and s90/— vl, the spectra are similar to those 
calculated using EoSQ. The pion V2{pt) is virtually insensitive to the change in freeze-out 
temperature, but the higher temperature leads to much larger pT^-differential f 2 for protons. 
This behavior has already been explained in Ref. [56|, where it was argued that the lower 
the temperature and larger the flow velocity, the smaller the V2{pt) at low values pt, and 
that the heavier the particle, the stronger this effect. Note that the three different lattice 
based parametrizations of EoS give almost identical results for V2 and spectra. This means 
that existing uncertainties in the EoS parametrization have negligible effect on the flow. 

Unfortunately it is not straightforward to calculate the spectra and f 2 using the various 
EoSs in the literature discussed earlier. One of the advantages of the Cooper-Frye procedure 
for freeze-out is that energy, momentum, particle number and entropy are conserved. But, 
they are conserved only if the equation of state is the same before and after the freeze- 
out [2^, i.e. that the fluid EoS is that of free particles and that the number of degrees of 
freedom in the fluid is the same than the number of hadrons and resonances which spectrum 
is calculated. This is not the case with the EoSs discussed here. When calculating the spectra 
we use the same set of resonances up to 2 GeV mass than what is included in our HRG. Laine 
and Schroder used resonances up to 2.5 GeV mass, but at, say, Tfo = 125 MeV freeze-out 



^ This temperature was found to reproduce the spectra when EoS Q is used 
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FIG. 9: The proton and pion spectra (left) and differential elliptic flow V2{pt) of protons and pions 
(middle and right) in 6 = 7fm Au+Au collisions for different EoSs. The results in the left and 
middle panels are calculated using the same freeze-out temperature Tfo = 125MeV for all the EoSs, 
whereas in the right panel it has been adjusted to produce similar p^-distributions. Tfo = 125MeV 
for EoSQ, and Tfo = UOMeV for EoSs s95p-vl, s95n-vl and s90/-vl. 



temperature the difference in energy and entropy is minuscule, about 0.05%. The situation 
with the other EoSs is more difficult. At the mentioned temperature, the Krakow EoS has 
4.5% smaller, EoSL 7% larger, and the HotQCD EoS 22% smaller energy and entropy than 
hadrons and resonances up to 2 GeV mass. One way to correct this discrepancy is of course 
to change the number of hadrons and resonances included in the spectra calculation. But, 
this approach would be tedious since the number of resonances needed to fit the densities in 
the EoS may depend on temperature. Also there is no telling whether there exists a set of 
resonances reproducing both densities and pressure at a given temperature, and the number 
of resonances could be surprisingly small. For example a hadron resonance gas consisting 
of only pseudo-scalar and vector meson nonets, and baryon octet and decuplet, still has 
10% larger energy density at T = 125 MeV than HotQCD EoS. Therefore we follow the 
approach espoused by Csernai 2^ and Bugaev 13] : We require that energy and momentum 



are conserved locally on the freeze-out surface, i.e. 

dcr^Tg^jj = dcr^Tpj^j,^i^[gg, (5.2) 

where T^^^^ is the energy-momentum tensor of the fiuid on the surface, and T^^^i^^^^icies 
energy-momentum tensor of the emitted particles. To conserve energy and momentum, 
we allow the temperature and flow velocity of fluid and particles differ, i.e. there is a 
discontinuity on the surface, and freeze-out is a shock-like phenomenon jSTJ. We have to 
admit the corrections due to this procedure are small and mostly affect the multiplicity, but 
consider obeying the conservation laws worth the extra effort. 

At flrst we use the same freeze-out energy density we used when comparing our 
parametrization to EoSQ, e = 0.065 GeV/fm^, for all EoSs. The corresponding temper- 
atures for each EoS are listed in Table IIIII In the left panel of Fig. [10] we show the pt 
distributions of pions and protons in 6 = 7 fm Au+Au collision. The differences in distri- 
butions are small, and the general behavior is what can be expected based on the stiffness 
of the EoSs and the flow anisotropy: The Krakow EoS is the stiffest, and leads thus to 
the flattest spectra. The EoS by Laine and Schroder leads to behavior almost identical to 
s95p-vl, and the HotQCD EoS and EoSL are the softest and have slightly steeper spectra 
than the other EoSs. For the pr-differential anisotropy of pions and protons the systematics 
is the same than seen for our parametrizations: when the freeze-out criterion is the same 
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efo 


^fluid 


-^particles 


efo 


^fluid 


-^particles 


s95p — vl 


0.065 


125 


125 


0.14 


140 


140 


HotQCD 


0.065 


129 


125 


0.11 


139 


135 


Laine 


0.065 


125 


125 


0.14 


140 


140 


EoSL 


0.065 


124 


126 


0.11 


134 


136 


Krakow 


0.065 


126 


125 


0.185 


146 


145 



TABLE III: The freeze-out energy densities (in GeV/fm'^) and corresponding temperatures (in 
MeV) for fluid and particles for each EoS used in the calculations shown in Fig. [10] 




FIG. 10: The proton and pion spectra (left) and differential elliptic flow V2{pt) of protons and 
pions (middle and right) in b = 7fm Au+Au collisions for s95p— vl (solid black line) and different 



EoSs in the literature 



24, 



4514471. |50| • The results in the left and middle panels are calculated using 
the same freeze-out energy density efo = 0.065 GeV/fm^ for all the EoSs, whereas in the right panel 
it has been adjusted to produce similar py-distributions, see text and Table Hill 



for all the EoSs, V2{pt) is basically independent of the EoS, as shown in the middle panel 
of Fig. [TUJ After the freeze-out criterion is adjusted to reproduce spectra obtained using 
EoSQ, the pion V2{pt) is independent of the EoS, but the proton anisotropy shows some 
sensitivity, see the right panel of Fig. [101 The differences between the EoSs are small and 
thus the differences in V2{pt) are small, but an ordering according to the stiffness of the 
EoS is visible: The Krakow EoS is hardest, and its proton V2{pt) is largest at small pt, 
whereas the HotQCD EoS and EoS L are softest and lead to lowest V2{pt) of protons at low 
Pt- After all the main results of this comparison are that the differences in the lattice EoS 
parametrization in the literature are small and not observable in the p^-- differential elliptic 
flow, and that energy conservation at freeze-out is not trivial if the EoS at freeze out is not 
that of free hadron resonance gas. 

Finally we want to compare the results of our calculations with data. Since all of our 
parametrizations lead to practically identical spectra and V2, we use only s95p-vl for sim- 
plicity, and compare the results to those obtained using EoS Q in Ref . [5| . First we fix all 
the parameters by requiring the reproduction of pion and net-proton [p — p) spectra in the 
0-5% most central Au+Au collisions at a/snn = 200 GeV energy. The resulting spectra 
are shown in the left panel of Fig. [Ill and the freeze-out temperatures are the same than 
mentioned before, Tfo = 125MeV for EoS Q, and Tfo = lAOMeV for EoSs s95p — vl. Since 
we assume chemical equilibrium we cannot reproduce both proton and anti-proton yields at 
such a low temperature. Our parametrization is for zero baryochemical potential, so we can- 



20 



Q. 
T3 



T3 




0.2 
0.16 
0.12 
0.08 
0.04 




Vp 




s95p-v1 




EoS Q 




STAR (open) 


/■^^^ 


PHENIX (filled)//' 










P 








Pt [GeV] . 



0.5 



1.5 



FIG. 11: Pion (vr^) and net-proton (p — p) spectra in 0-5% most central (left), and pion and 
antiproton p^-difFerential elliptic flow vi(j>t) in minimum bias (right) Au-|-Au collisions at ^snn = 
200 GeV compared with hydrodynamic calculations using two different EoSs and assuming chemical 
equilibrium. The data was taken by the PHENIX 58 1 and the STAR [s^ collaborations. 



not calculate net-protons either, but we approximate them by having a finite baryon density 
in the calculation, and converting this density into a finite chemical potential at freeze-out 
by using a HRG EoS which allows a finite net baryon density. The right panel of Fig. [TT] 
shows the pr-differential elliptic flow of pions and anti-protons in minimum bias Au-|-Au 
collisions at the same energy. As earlier, the pion V2{pt) is very similar for both EoSs, but 
the anti-proton V2{pt) is quite different. In fact for the realistic EoS the anti-proton V2 is 
largely overpredicted, while we have a reasonable agreement with the data when the bag 
model EoSQ is used. This is very similar to the finding of Ref. (sf. Also the Krakow group 
seem to overpredict the proton V2 46| indicating that this could be a general feature of ideal 
hydrodynamic models using more realistic EoS. The first results indicate that dissipative 



effects can at least reduce this problem [60 



This analysis has been repeated assuming partial chemical equilibrium in the hadronic 
phase. The details of this analysis are given in the appendix Our main finding did not 
change. Using different EoS for the same initial conditions and kinetic freeze-out temperature 
has little effect on w 2 but leads to significant differences in the spectra. If the spectra are fitted 
to reproduce the experimental data by adjusting the initial conditions, the pr-differential 
elliptic flow of anti-protons becomes too large for lattice based EoS, but agrees reasonably 
well with EoS Q. On the other hand, the p^-differential elliptic flow of pions is clearly too 
large for both EoSs requiring substantial dissipation to reduce it to fit the data. 



VI. CONCLUSION 

In this paper we addressed the question to what extent the Hadron Resonance Gas (HRG) 
model can describe the thermodynamic quantities calculated on the lattice. As discussed 
above this question is very important for implementing a consistent freeze-out prescription 
or a consistent switch to transport in hydrodynamic or hybrid models, respectively. We 
have found that lattice data strongly disagree with the HRG model in the low temperature 
regime. The reason for this disagreement has been identified with large cutoff effects in the 
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lattice calculations. We also showed that taking into account the discretization effects in 
the hadron spectrum in the HRG model leads to a good agreement with the lattice data. 
In fact, we find that for some quantities the HRG model works well to unexpectedly high 
temperatures. Based on this observation we constructed several parametrizations of the 
equation of state which interpolate between the lattice data at high temperature and the 
resonance gas in the low temperature region. The central quantity in this analysis was the 
trace anomaly since it is directly calculated on the lattice and the differences in the proposed 
parametrizations are found in the temperature region where the trace anomaly reaches its 
maximal value. 

We studied the hydrodynamic evolution using three parametrizations of the EoS that 
interpolate between HRG EoS and the lattice data and compared the results with the corre- 
sponding ones obtained using an EoS with a first order phase transition, the so-called EoS Q, 
as well as several other parametrizations of the EoS used in the literature. We have analyzed 
the flow in terms of momentum space anisotropy e^, p7-differential elliptic flow V2{pt) and 
proton and pion spectra. The three parametrizations of the EoS proposed in this paper as 
well the the parametrization by Laine and Schroder [i^] gave very similar results for all of 
the above quantities. The effect of using different EoS parametrizations is the most visible 
in ep. The difference in the results obtained with EoS Q and other parametrizations is espe- 
cially large. Quite surprisingly V2{pt) is not sensitive to the choice of the EoS if the same 
freeze-out temperature is used. The particle spectra on the other hand are sensitive to the 
EoS. However, the change in the EoS can be compensated by change of the freeze-out tem- 
perature. If the freeze-out temperature is adjusted to reproduce the particle spectra we see 
large differences in the proton V2{pt) for EoS Q and other EoS parametrizations. However, 
for all the other parametrizations considered here, the proton V2{pt) is quite similar. 

The work presented in this paper should be extended in number of different ways. First 
we should extend the comparison of lattice QCD results with modified HRG to other fluctu- 
ations, including electric charge fluctuations as well as fourth and higher order fluctuations 
of baryon number, strangeness and electric charge. However, since these quantities were 
studied in detail only for the p4 action, in the analysis additional assumptions about the 
cutoff effects in the hadron spectrum have to be made. Furthermore, it would be interesting 
to study the quark mass effects in the HRG model, especially since recent lattice calcula- 



tions of the EoS extend to physical values of the light quark masses [61] ■ We also should 
consider the effect of finite baryon potential on the EoS. This will become important for 
application of hydrodynamic models to heavy ion collisions at lower energies, especially to 
the proposed RHIC energy scan. Finally, it will be interesting to study the effect of the EoS 
in the framework of viscous hydrodynamics. We plan to address these issues in forthcoming 
publications. 
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1.51418 
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n 
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TABLE IV: The values of the parameters appearing in Eq. ()Aip for different hadrons. 

for comparison. P.P. thanks the members of the MILC collaborations, especially Carleton 
DeTar, Steve Gottlieb, Urs Heller and Bob Sugar for correspondence and for providing their 
numerical results, including the unpublished data of Ref. |37|. P.P. is also grateful to Zoltan 
Fodor and Sandor Katz for correspondence and for providing the pion masses for the stout 
action. 



Appendix A: Hadron masses on the lattice 

In this appendix we are going to discuss the cutoff and quark mass dependence of hadron 
masses calculated on the lattice and give the parameters entering Eqs. fl2.8p - fl2.12p . We 
have fitted the quark (pion) mass and lattice spacing dependence of the p, K*, (p, N and 
Q masses obtained in Refs. [25-2^ by a simple Ansatz 



aiinm^f bix / / n2 /atn 

rim = rirriQ H — ; h ; — , x = (a ri) (All 

1 + a2X 1 + 02X 



The values of the fit parameters mo, ai, 02, hi and 62 are given in Table HVl In Figs. [121 
and [H] we show the above parametrization against the available lattice data. Here we note 
that the lattice data for the VL mass have been corrected to take into account that the physical 
strange quark mass is slightly smaller that the one used in lattice simulations. It turns out 
that Eq. flAip reproduces the experimental values of the hadron masses in continuum limit at 
the physical point rim^ = 0.226. This justifies the use of Eqs. f l2.8p - fl2.12p for the evaluation 
of the hadron masses in the HRG model. 

As discussed in the main text due to lack of detailed lattice studies we used Eq. (12. 8p 
and Eqs. fl2.10p - fl2.12p to evaluate the mass of the A resonance as well as single and double 
strange baryons with values of the parameters in Table [IV] corresponding to the nucleon. In 
Figure [15] we compare our estimates of the A, A and S baryon masses shown as lines with 



available lattice data from the MILC collaboration 1271. 28 



Appendix B: Fitting procedure 

Here we describe in detail how we fit the trace anomaly of lattice and hadron resonance 
gas. The two first terms of the inverse polynomial Ansatz (Eq. 14.20 

(e - 3P)/r^ = d2/T^ + d^/T^ + ci/T"! + Ci/T''^ (Bl) 

appear to provide good fit of the lattice data at high temperatures, T > 250 MeV ^23]. We 
want to join this parametrization to the trace anomaly of hadron resonance gas and require 
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FIG. 12: The p-meson (right) masses calculated on the lattice using asqtad action 271. l28l. |37I| and 
compared with Eq. ()Aip (lines). 
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FIG. 13: The K* (left) and 0-meson (right) masses calculated on the lattice using asqtad action 
27|, |28|] and compared with Eq. (|Aip (lines). 



that the trace anomaly and its first and second derivative with respect to temperature are 
continuous where joined. Thus, we need one additional term with negative coefficient ci 
and exponent ni > 4 to produce a peak around T ^ 200 MeV, and another with positive 
coefficient C2 and exponent n2 > rii to make the second derivative continuous. We calculate 
the trace anomaly of hadron resonance gas using all the resonances up to 2 GeV mass^ 
in the summary of the 2004 edition of the Review of Particle Physics [63]. We note that 



The list of included resonances and their properties can be found at [62[ . 
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FIG. 14: The nucleon mass 271. l28l. l37l| (left) and the O baryon mass 25|] (right) calculated with 
asqtad action and compared with our parametrization. 
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FIG. 15: The baryon masses calculated for asqtad action [27|, l28|, l37| and compared with our 
parametrization. 



including all the resonances up to 2. 5 GeV instead of 2 GeV mass gives result for the trace 
anomaly, which is only 2% larger at T = ITOMeV and 3% larger at T = ISOMeV. This 
change is definitely smaller than the expected discrepancies between HRG and lattice at 
these temperatures. 

In the ansatz we have seven unknown parameters: the coefficients ^2, c?4, C\ and C2, expo- 
nents n\ and n^^ and the switching temperature Tq. We have four constraints, the continu- 
ity of the trace anomaly and its derivatives at Tq, and the requirement s(T = 800 MeV) = 
0.95-ssB or s(T = 800 MeV) = 0.90- ssb- The requirement of the continuity of the derivatives 
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gives two equations to fix the parameters c\ and C2- 



n2(?T.2 - rii) ^2(^2-^1) 
ni(?7,i-r;,2) ni(ni-n2) 



where 



d e-W 



dT T4 



G'i(T) and 



rf2 e-3P 



From Eg. (14. II) and Ts = e + P we obtain 



= G,{T). 

HG 



± = d ^ A _ ^ V ^ : ^ M , + ^ r A : : 4P(To) 

Since entropy at T = 800 MeV is fixed, we can use the above equation to constrain d^, 

di = d4{d2, ni, n2, Tq). (B3) 

We can thus express the parameters ci, C2 and ^4 in terms of ^2, ni, n2 and Tq, and use the 
continuity of the trace anomaly to fix Tq. We get an equation 



e-3P 



2^4 



/^^_c^2 , d4{d2,ni,n2,TQ) ci{d2,ni,n2,To) 02(^2, ni, n2, Tq) 
UoJ-;^H Tf^A ^ ^ ^^^^i ' y^^) 

HG -'O -'O -'O -'O 



which can be evaluated numerically to obtain Tq. This procedure leaves us with three un- 
knowns d2, rti and n2, which are chosen to fit the lattice data. However, such a fitting 
procedure would be highly nonlinear. We simplify the problem by requiring that the expo- 
nents are integers, and use brute force: We make a single parameter ((^2) fit with all the 
integer values 4 < ni < 31, and rii < n2 < 43, and choose the values ni and n2 which lead 
to the smallest x^- Alternatively we can fix Tq to a prescribed value, use Eq. ( 1B4I) to fix the 
value of d2, and use only rii and n2 to perform the fit. 

In our fit we use the lattice data for T > 250 MeV obtained with p4 action on A^^ = 8 
lattices as they extend to sufficiently high temperature In addition, we include = 6 
p4 data for T > SOOMeV joj. The fits in general do not reproduce the lattice data in peak 
region (190MeV < T < 250MeV). On the other hand the height of the peak in the trace 
anomaly may be affected by discretization effects. This can be seen as a difference between 
the Nr = 6 and iV^ = 8 results. Assuming that discretization effects in the peak height go 
like we can estimate the trace anomaly at T = 206MeV to be 5.7 ± 0.15. We can use 

this value as an additional data point in our fits. We label different parametrization of the 
trace anomaly obtained using different constraints on the entropy density at T = 800MeV 
and its height in the temperature region IQOMeV < T < 250MeV as s95p— vl, s95n— vl, 
s90/— vl etc. The first three characters stand for the constrain on the entropy density (90% 
or 95% of the ideal gas value). The fourth character stands for additional constraints on 
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the trace anomaly in the peak region. Namely, "n" stands for a fit with no constraints: 
using data for T > 250 MeV, and Tq as a free parameter in the fit. "p" means having an 
additional data point of 5.7 ± 0.15 at T = 206 MeV to constrain the peak, and "f" stands 
for a fixed value Tq = 170 MeV in the fit. Finally "vl" is the version number of the current 
parametrization. The value of the parameters ^2, (^4, Ci, C2, ni, n2 and Tq for different fits 
are given in Table [ITl This procedure was designed for numerical applications when the trace 
anomaly is numerically evaluated using Eq. fl2.ip and laws of thermodynamics. For practical 
purposes we also provide a parametrized version of the trace anomaly of the hadronic part 
of our EoS. We choose a polynomial 

e — 3P 

= aiT'i + asT'^ + agT'^ + a^T^^ (B5) 

and fit it to the trace anomaly of the hadron resonance gas evaluated in the temperature 
interval 70 < T/MeV < 190 with 1 MeV steps assuming that each point has equal "er- 
ror" . The limits have entirely utilitarian origin: in hydrodynamical applications the system 
decouples well above 70 MeV temperature and only a rough approximation of the EoS, 
P = P{e), is needed at lower temperatures. On the other hand we expect to switch to the 
lattice parametrization below 190 MeV, and the HRG EoS above that is not needed either. 
We fix the exponents in Eq.( ]B5p again using brute force. We require them to be integers, 
go through all the combinations < /i < ^2 < ^3 < ^4 < 10, fit the parameters ai, a2, as, 
to the HRG trace anomaly evaluated with 1 MeV intervals, and choose the values /i, I2, I3 
and U which minimize the x^- We end up with li = 1, I2 = 3, I3 = 4, /4 = 10, and ai = 4.654 
GeV-\ a2 = -879 GeV'^, = 8081 GeV-^ 04 = -7039000 GeV"^" . To obtain the EoS, 
one also needs the pressure at the lower limit of the integration (see Eq. fl4.1l) ) Tiow = 0.07 
GeV: P(Tio„)/Tj^^ = 0.1661. Our EoSs are also available in a tabulated form at |62| . 



Appendix C: Spectra and elliptic flow for partial chemical equilibrium 

In this appendix we discuss the elliptic flow and the spectra of protons and pions and their 



sensitivity on EoS when partial chemical equilibrium |44l . |64| . |65| is assumed to reproduce the 
observed particle yields. The EoS for the system in partial chemical equilibrium is available 
in tabulated form [62] . We again calculate the flow in Au+Au collision at ^snn = 200 GeV 
with impact parameter 6 = 7 fm. First we used the same initial condition and the same 
freeze-out temperature for all EoSs, namely Tchem = 150MeV for the chemical freeze-out 
temperature and Tkin = 120MeV for the kinetic freeze-out temperature. The initial time 
for the hydrodynamic evolution was chosen to be r = 0.6fm and the initial entropy density 
was chosen as for the case of chemical equilibrium. The results are shown in Figure [161 As 
one can see the elliptic flow is not very sensitive to the choice of EoS if everything else is 
kept unchanged, and all three parametrizations used in the analysis give almost the same 
result. The spectra are much more sensitive to the EoS but again there is no difference for 
the different lattice parametrizations. 

For a proper discussion on the sensitivity of elliptic flow on the EoS, one has to again 
re-tune the calculation to reproduce the experimental results on the pj^-spectra. We 
obtained the best results by keeping the chemical and kinetic freeze-out temperatures, 
^chem = 150MeV and Tkin = 120MeV unchanged, and tuned the initial conditions. For 
EoS Q we used tq = 0.2fm and initial entropy density which scales with the number of 



binary collisions (see Ref. 4J]). For the s95p-vl parametrization of the EoS we used two 
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FIG. 16: The p^-differential elliptic flow V2{pt) of protons and pions (left) and proton and pion 
spectra (right) for different EoSs in 6 = 7 fm Au+Au collision when chemical freeze-out takes place 
at Tchem = 150 MeV and kinetic at Tkin = 120 MeV. 



different initial conditions. One with tq = 0.8fm and initial entropy density proportional 
to the number of binary collisions, and a second one with tq = 0.2fm and initial entropy 
density scaling with a combination of binary collisions and number of participants. The 
corresponding results are shown in Figure [T71 We see again that the lattice based EoS give 
larger pr-differential elliptic flow for the protons than EoS Q for both initial conditions. In 
this case, however, the EoS Q does not do a good job of describing the data either, in agree- 
ment with previous findings |4i,|65|. Especially the pion p-r-differential V2 is too large for 



all EoSs, and there is clearly room for significant dissipation to reduce the anisotropy. It is 
also worth to notice that the uncertainty related to the initial state is at least as large as 
the effect of the EoS on the proton V2{pt) and better theoretical constraints to the initial 
state are needed. 
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